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ABSTRACT 

It is shown that the least squares collocation approach to estimating geodetic parameters 
is identical to conventional minimum variance estimation. Hence the least squares collocation 
estimator can be derived either by minimizing the usual least squares quadratic loss function 
or by computing a conditional expectation by means of the regression equation. 

When a deterministic functional relationship between the data and the parameters to be 
estimated is available, one can implement a least squares solution using the functional relation 
to obtain an equation of condition. It is proved the solution so obtained is identical to what 
is obtained through least squares collocation. The implications of this equivalance for the 
estimation of mean gravity anomalies are discussed. 
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ON LEAST SQUARES COLLOCATION 


INTRODUCTION 

A characteristic of geodetic research is that numerous data types are available for esti- 
mating parameters of interest. The problems of combining heterogeneous geodetic data types 
to provide consistent estimates has lead some researchers to the belief that conventional least 
squares methods are inadequate. An alternative approach to geodetic data reduction prob- 
lems called least squares collocation has been suggested by Moritz [ 1 ] . Some authors have 
claimed that least squares collocation is a more general and more powerful parameter estima- 
tion procedure than the classical least squares method [ 1, 2, 3, 4, 5] . It has also been asserted 
that least squares collocation is the only parameter estimation method which permits the 
simultaneous and optimal processing of heterogeneous data types [6, 7] . Other authors have 
disputed these claims [8, 9] . 

This note is an effort to settle what has become a confusing and contentious issue. It 
will be demonstrated that least squares collocation is an estimator of a type which is well 
known in conventional estimation theory. The presentation is elementary in content and 
should be intelligible to anyone familiar with the rudiments of probability theory. 
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SOME PROPERTIES' OEiMINIMUM VARIANCE ESTIMATORS 'f 
Let X be a finite dimensional.-vector of parameters to be estimated. Sinc,e:?the’param'; 
eters are not perfectly known it is legitimate to view X as a random vector. Also there is- no ■ 
loss, in generality in assuming the zero vector to be the expectation of X; Let the covariance 
matrix of X be known. Thus: 

E(XXT) = C (1) 

where C is positive definite. Assume the existence of a finite vector Y which defines a state 
which is directly observable. Hence Y is a random vector which is sampled by a measuring 
process. . ' 

.. Lacking data, the minimum variance estimate of X is the zero vector. Butdntuitively it 
> is clear that if random vectors Y and X are correlated and if a realization Y' of Y is, available, 

.A 

it should be possible to obtain an improved estimate X of X. Several, criteria are available; 
Two of the most commonly used are 

Criterion A - choose X as that vector which minimizes the conventional least 
squares quadratic form. 

A 

Criterion B - choose X as the expectation vector of the conditional distribution 
of X given a realization Y' of Y. 

It will be shown that the application of either criterion leads: to the same estimator. 

A 

To' obtain the improved estimate X,,it:is necessary to precisely define the correlation 
between Y and X. This is commonly donetin two ways which we will describe as a model 1 
and model 2. In model 1 the correlation is, described by a linear stochastic equation, 

Y = SX + i^,?E(i^) = Oi;E(wT) = Q (2) 

In model 2 the correlation is described in. terms of a^eross covariance matrix. 
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In fact, models 1 and 2 are alternative and equivalent ways of describing the second order 
statistical properties of the joint distribution of Y and X. Model 1 can be transformed to 
model 2 by defining the symbols A and B on the right side of equation 3 as 

A - S C S'^' + Q B = SC (4) 

Conversely, model 2 can be converted to model 1 as described in equation 2 with * 

S = BC-i Q = A-BC-’B^ (5) 

Arbitrarily, we choose model 1 as a description of the necessary correlation. The applica- 
tion of criterion A imphes the minimization of the quadratic loss function 

L(X) = (Y' - SX)'^ Q*i (Y' - SX) + X'^ C*’ X (6) 

where Y' is a realization of Y. The solution to the minimization problem is 

X = (S'f Q-’ S + C*’)"’ (S^ Q-'' Y') (7) 

Equation 7 is known to represent a minimum variance estimator [ 10] . 

To apply criterion B, transform model 1 to model 2 by means of equation 4. The 
well known regression equation [10] can then be employed on the right side of equation 
3. The result is 

X = E(XlY = Y') = B^ A-’ Y' (8) 

Again using equation 4, equation 8 can be transformed into 

X = C S^(S C S'^-l- Q)-’ Y' (9) 

The Shure matrix identity can be used to translate equation 9 into the alternative form: 

X = (S'^ Q-’ S + C'^ )■’ Q-'' Y' (10) 

Equation 10 is identical to equation 7. We have established the following 

Theorem 1 - Assume that Y and X are correlated random vectors and that a 
realization Y' of Y is available. The correlation may be defined either in terms 
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of a linear stochastic equation given by equation 2 or a cross covariance .matrix 
given by equation 3 . In each case the minimization of the least squares ^quadratic 
loss function of equation 6 and the computation of the expected value of the 
conditional distribution of X given Y' by means of the regression equation yield 
the same minimum variance estimator. 
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LEAST SQUARES COLLOCATION 

Let Y' be a set of geodetic observations. The problem is to obtain from such an obser- 
vation set an estimate of a set of geodetic parameters X. The starting point of the least squares 
collocation approach to the problem is the assumption that one has full knowledge of the 
second order statistics of the anomalous potential. Let P (xi) and P (X 2 ) be the anomalous 
potentials at points xi and xj on or outside the reference geoid. We assume the possession of 
a function K (xi, X 2 ) such that 

E(P(Xi)P(X2))=K(Xi,X2) (11) 

Equation 1 1 defines the so-called covariance function. Let £ be the countibly infinite set of 
deviations of the spherical harmonic coefficients of the Earth’s potential field from reference 
values. A convenient way to define a covariance function is to specify the second order 
statistics of £. Hence define 

E(jejC’’) = T (12) 


The matrix T uniquely defines a covariance function. Algorithms for determining the right 
side of equation 1 1 given the right side of equation 12 may be found in Moritz [ 1 ] or 
Tscherning and Rapp [11]. Conversely, a given covariance function uniquely defines a 
covariance matrix T [12] . Hence there is no loss in generality in assuming that the co- 
variance function for the least squares collocation procedure is given in terms of a ma- 
trix T as defined by equation 12. Let Y be the ideal observation state of which Y' is a reali- 
zation. Since both Y and X are geodetic entities they are functions of the set of spherical 
harmonic coefficients of the Earth’s potential. First order Taylor series expansions of these 
functions about a set of reference spherical harmonics will yield linear matrix equations 


a, Y = /iT 

b, X = /2 jC 


(13) 


where reference values of Y and X are assumed equal to the zero vector. In equations 1 3 and 
in subsequent equations, whenever the matrix symbolism implies countible infinite summa- 
tion it is the limiting value which is intended. Alternatively, the reader can assume that the 
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vector £ of deviations of spherical harmonic coefficients from referetice valuestis truncated 
at a sufficiently high degree that errors in representation in equations 13 are negligible. 
Equations 1 2 and 1 3 yield 

a, E(YYT) = A = /iT/y 

b, E(XXT) = C = /2T/^ (14) 

c, E(YXT) = B-/iT/T 

The actual observations Y' are corrupted by errors in the measuring system. Hence 

Y' = GZ + Y+i; (15) 

where 

a, 
b, 
c, 
d, 


E(j;) = 0, 

E(wT)^q 

E(z) = 0, 
E(zzT) = P, 




(16) 


The vector Z is interpreted as a set of parameters which determines the systematic part of the 
errors in the measuring system. Define an augmented parameter set as 

(17) 

Equations 14, 15 and 16 define the correlation between random vectors Y' and S. In the 
previous section it was demonstrated that given the correlation between two random vectors 
and given a realization of one of the vectors, it was possible to construct a minimum variance 
estimator for the other vector. This estimator can be obtained either by using the regression 
equation to compute a conditional mean or by minimizing a conventional least squares quad- 
ratic loss function. Arbitrarily we will obtain the minimum variance estimate -for S by comput- 
ing the conditional mean of S given a realization of Y'. The covariance matrix for the joint 
distribution of S and Y' is 
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E 


(18) 


where 


Y' 

S 


Y' 

S 



A + G P G'^ + Q 
PT 


F 

D 


a, F = [b,GP] 


b, 



0 

P 


(19) 


Let Sc be the conditional distribution of S given a realization of Y'. By assuming either that 
the random vectors are normally distributed or that the expectation of Sc is a linear function 
of the measurements, we can resort to the regression equations for the mean and co- 
variance of Sj. as 

E(Sc) = S = F'T(A + GPGT + Q)-i Y' (20) 


E ^(Sc - S) (Sc - S)'^) = D - pT (A + G P G"^ + Q)-’ F (21) 

Equations 17, 19, and 20 permit us to separately write the conditional expectations of X and 
Z as 


a, X = B"^ (A + G P gT -I- Q)-’ Y' 

b, Z = P G"^ (A + G P G'f + Q)*^ Y' 


( 22 ) 


A straightforward application of the Shure matrix identity converts equations 22 into 


a, X- B’f(A + Q)-^Y' -G Z) 
b, Z = (G(A + Q)-’ G"^ + P-’)*'' G"^ (A + Q)'’ Y' 


(23) 


Equations 23 represent the least squares collocation solution for geodetic paramv'ters X and 
measuring system parameters Z given a realization of observation vector Y'. We have proved 
the following. 

Theorem 2 : The least squares collocation splution for geodetic parameters X and 
measuring system parameters Z given a realization of an observation vector Y' is 
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identical to the conventionatoihimum variance solution. Hence the eolfecation 
solution can be obtained either by determining the conditional expectations of X 
and Z by means of the regression equation or by minimizing the usuaf least squares 
quadratic loss function. 

Tapley [8] provides a somewhat longer proof of the equivalence of least squares collo- 
cation to the conventional minimum variance estimator. Tapley ’s proof is interesting because 
it relies entirely on elementary matrix operations. 



APPLICATIONS 


The previous sections show that the techniques of conventional minimum variance esti- 
mation have considerable power and generality. With the appropriate use of these techniques, 
one can obtain a minimum variance estimate of any set of parameters which are functions 
of the anomalous potential given a realization of any observation set which is also a function 
of the anomalous potential. In many cases it is possible to augment the parameter set in 
question in such a way that the laws of Mathematical Geodesy provide a deterministic func- 
tional relationship between the augmented parameter set and the ideal or noiseless represen- 
tation of the data set. For this situation two different estimation procedures are available: 


Estimation Procedure 1) Use the postulated covariance matrix for the anomalous spherical 
harmonic coefficients to construct the covariance matrix for the joint distribution of the 
data set and the parameter set to be estimated. The regression equation can then be used to 
compute the conditional mean of the parameter set. (least squares collocation) 


Estimation Procedure 2) Use the postulated functional relationship between the augmen- 
ted parameter set and the data set to construct an equation of condition for a conven- 
tional least squares with a priori estimation of the parameter set. 

It will be shown that the two procedures are equivalent. Again let T represent the co- 
variance matrix for the anomalous spherical harmonic coefficients £. Let Y be the ideal or 
noiseless representation of the data and let Xibe the parameter set to be estimated. Assume 
that Xi is part of a larger parameter set. Hence 



(24) 


The vectors Y and S are functions of the anomalous potential. First order Taylo i- series ex- 
pansions of the appropriate functions yield linear matrix equations 
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( 25 ) 


a, Y = U£ 

b, Xi=/2,i£ 


C,. X2 - /2,2 £ 


Equations 25b and 25c can be rewritten as 


S = /2X 


where 


f2- 


/2,1 

/2,2 


(26) 


(27) 


For simplicity assume that there are no systematic errors in the measuring system. Then'the 
observation equation is 

Y' = Y + v, E(j^) - 0, E(wT) = q (28) 

Equations 25 and 28 permit us to write 

a, E(Y'Y^) = /,T/'[+Q 

, T T (29) 

b, E(Y'x}) = /iT/Ji 

The least squares collocation solution for Xi is 

Xi - /2,i T,/y (/i T /[ + Q)*’ Y' (30) 

To implement estimation procedure 2, assume that the laws of Mathematical Geodesy 
provide a deterministic functional relationship between Y and S. A first order Taylor series 
expansion of the function will yield 


Y = /3 S 

Notice that from equations 25a, 26 and 31 we have 

/l = /3/2 


(31) 


(32) 
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Given equations 28 and 3 1 , the usual least squares solution for S is given by 

S = {fs Q-'fl + (/2 T fl )' y fl Q-’ Y' (33) 

Transform equation 33 to the equivalent form 

s = /2 T /T /T (/3 /2 T /J/| + Q)-' Y (34) 

Withe the aid of equation 32, equation 34 can be rewritten as 

S = /2T/T(/,T/j + Q)-iY' (35) 

Hence 

Xi = /2.1 T/T(/,t/T + Q)-1y' (36) 

Equation 36 is identical to equation 30 and this proves the equivalence. 

The results of this section have important consequences for the estimation of mean 
gravity anomalies. Stokes’ formula provides a representation of the anomalous potential on 
or outside of a reference geoid as [ 1 3] 

U(r,0,X) = ||:||s(r,.//)6gda (37) 

where r, 0, and X are the spherical coordinates of the computation point, R is the value of the 
Earth’s radius, S(r,0) is the Stokes function with 0 the spherical distance between the inte- 
gration point and the projection of the computation point on the reference geoid. The 
symbol 5g is the point gravity anomaly referenced to the nominal field and measured on the 
reference geoid. The integration is over the entire geoid. The discrete approximation to 
Stokes’ formula is 

U(r, 0, ^ 2 S(r,0i) dOi (38) 

1 

where 5gj is interpreted as a mean gravity anomaly averaged over a non-zero surface area 
da,. The summation of equation 38 is finite and encompasses the entire reference geoid. 

Let Y be a vector which is determined by the anomalous potential field. Assuming the 
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validity of the discrete form offStokes’ formulaj equation 38 provides a relationship be- 
tween Y and a globally distributed set 6g of mean gravity anomalies which;; after suitable, 
linearization can be written as 

Y = /% (39) 

Equation 39 can be used as an equation of condition for a least squares with a priori 
estimate of a global set of mean gravity anomalies. On any subset of the globally dis- 
tributed mean gravity anomalies the solution so obtained will agree with the least squares 
collocation solution. 



COMMENTS 


The least squares collocation algorithm can be exhibited as a conventional minimum 
variance estimator. Hence, the algorithm can be derived either as an application of the 
regression equation or by minimizing the usual least squares quadratic form. In some 
cases a geodetic parameter set to be estimated can be augmented in such a way that the 

I 

laws of Mathematical Geodesy provide a deterministic relation between the augmented 
parameter set and the available data. For this case the deterministic relation can be used 
to obtain an equation of condition for a conventional least squares with a priori estimate. 
The solution so obtained must agree with the least squares collocation solution. 

For estimating mean gravity anomalies, both least squares collocation and the con- 
ventional least squares approach utilizing Stokes’ formula are applicable. Each procedure 
must employ a certain approximation. With the conventional least squares approach an 
integral representation (equation 37) is replaced by a finite sum (equation 38). With least 
squares collocation, covariance and cross covariance representations for point gravity 
anomalies must be averaged to obtain similar representations for mean gravity anomalies. 

In each case the approximations can be performed so that the errors of representation 
are less than any preassigned value. The results of this paper show that if the two esti- 
mation procedures are implemented in such a way that corresponding errors of representa- 
tion are negligible, resulting estimates of mean gravity anomalies will be equal. Hence, 
the choice between estimation procedures should be mode on the basis of computational 
convenience. 

A disadvantage of the conventional least squares approach to estimating mean gravity 
anomalies which relies on a discrete form of Stokes’ formula is that its rigorous implemen- 
tation implies the simultaneous estimation of a global set of anomalies. With the least 
squares collocation approach it is convenient to estimate anomalies on a one by one ba- 
sis. However, it can be shown that for many data types and with proper estin.ation strat- 
egies [14, 15, 16, 17] , it is possible to estimate local blocks of mean gravity amimalies 
without serious aliasing. 

A serious computational problem associated with least squares collocation is that its 
implementation implies the inversion of a matrix whose dimension is the size of the data 
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se’t. The conventional least squaref ip^^^^ implies the inversion of a* matrix Mii^se dimen- 
sion is the size of the parameter set to be estimated. Hence, when large and dense data 
distributions are available for estimating mean gravity anomalies a conventibhai least squares 
technique utilizing Stokes’ formula is a more lo^cal choice for an estimation procedure. 
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